%% This code runs Stage III of the method.
% assuming constant cost
close all; clear;
clc;

% define saved files from running this code.
estFile_save = 'est2_output1.mat';
fig1_save = 'fig3_supply_mean.png';

%% load data and estimates from Stage I & II
load 'ATE_data_age10.mat'
load 'est1_output_2015.mat' % output year: 1950-2006

Tend = 2015;
ttmax = length(E);
H = (1e+3)*ATE_catch(year<=Tend); % harvest (kg)
Price = R16price(year<=Tend); % price ($/kg)

%% Estimate coefficients in the dynamic Open Access difference eq
fprintf('DeltaE = gamma*(PH) + (gamma*c)(-E)\n');
% c(t) = c, costant unit cost of effort

% generate variables for regression y = X*b
y = E(2:end)-E(1:end-1); % DeltaE
z = Price(1:end-1).*H(1:end-1); % Price*Catch
x0 = -E(1:end-1);

X = [z x0];
fprintf('y = DeltaE, x1 = PH, x2 = (-E)\n');
result = fitlm(X,y,'Intercept',false);
disp(result)
[n, p] = size(X);
AIC = n*log(result.SSE/n)+2*p;

b1 = result.Coefficients.Estimate; % point estimates
se1 = result.Coefficients.SE; % std. err
gamma = b1(1);

%% save residuals (for appendix Fig E.1)
% resid1 = y - result.Fitted;
% save('est2_output1_residual_2015.mat','resid1')
%% Monte Carlo 
M1 = 1000; % number of random draws
M2 = 500; % number of MC iterations

rng(54321) % for replication

% random sampling of gamma from Gamma dist.
pd1 = makedist('Gamma','a',(b1(1)/se1(1))^2, 'b',se1(1)^2/b1(1));
beta1_mc = random(pd1,M1,M2); 
% random sampling of gamma*c from Gamma dist.
pd2 = makedist('Gamma','a',(b1(2)/se1(2))^2, 'b',se1(2)^2/b1(2));
beta2_mc = random(pd2,M1,M2);
% sample mean of cost parameter
c_sample = mean(beta2_mc./beta1_mc);

c = mean(c_sample); se_c = std(c_sample);
%% save cost and gamma estimates
save(estFile_save,'gamma','c','se_c')

%% Simulate backward-bending supply using c
cost = c; Rscale = R1_scale;

[H_sim,P_sim,Hmax,backP] = sim_supply(Rscale,q0,Theta,cost);

fig1 = figure('Position',[0 0 700 500]);
hold on
plot((1e-3)*H_sim,P_sim,'k','LineWidth',2);
line([0 (1e-3)*Hmax],[backP backP],'Color','k','LineStyle',':','LineWidth',2)
xlabel 'Catch Quantity (metric tons)'
ylabel 'Price ($/kg)'
ylim([0 100])
xlim([0 Inf])
%gtext(sprintf('$%.2f/kg',backP),'FontSize',12)

saveas(fig1,fig1_save)